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Bessel functions with pure imaginary index (order) play an important role in corpuscular optics 
where they govern the dynamics of charged particles in isotrajectory quadrupoles. Recently they 
were found to be of great importance in semiconductor material characterization as they are man- 
ifested in the strain state of crystalline material. A new algorithm which can be used for the 
computation of the normal and modified Bessel functions with pure imaginary index is proposed. 
The developed algorithm is very fast to compute and for small arguments converges after a few 
' iterations. 
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Bessel functions occur in many branches of mathematical physics as solutions of differential equations when boundary 
' conditions such as the Dirichlet or Neumann are imposed on various space domains. The Bessel functions with index 
| (order) v can be represented as the solutions of the following differential equation [l[ : 
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I. INTRODUCTION 



x 2 y" + xy' + (x 2 ~v 2 )y = 0, (1) 

where y' and y" are the first and second derivatives (respectively) with respect to x while v is a complex constant. 
The solutions of eqn.{T]) can be expressed as an absolute converging series that is defined in the entire complex plane: 

!2: J ^ = ®"t(-v nT{v l +n+1 ". m 

, In a great variety of applications such as diffraction, electrical induction, etc, only Bessel functions of the zeroth 
■ and first orders, namely, Jq(x) and J\{x) occur while in other physical applications (such as solutions of Kepler's 

t-H ' equation), the entire indices n are used. Bessel functions with complex indices were considered to have limited areas 

Q\ \ of applicability in natural and applied sciences till recently. 

f^) . The normal and modified Bessel function with pure imaginary index was shown to be the solution governing the 
motion of charged particles in isotrajectory quadrupoles 0. Recently these functions have been found to be present as 
solutions of the Lame's equation which characterizes the displacement (and strain) field distribution in semiconductor 
nanostructures Q- 

In order to compute the solutions of equation ([T]) for the case of purely imaginary indices if, it is required to obtain 
the sum of a series with each coefficient being a complex number and to compute the complex value function T{iv) 
which on its own is a very laborious operation. In reality, for the case of a purely imaginary index iv and for a natural 
number n we have: 

n 

T{iv + n + 1) = T(w) Y[ {w + m) (3) 

m=0 



and thereafter Bessel function @ may be re-expressed in the form 
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where v and x are real and complex numbers respectively. The solution of JTJ) expressed in the terms of (|4]) represents 
a "dimcult-to-computc" complex solution of a differential equation with the real coefficients. To the best of our 
knowledge there is no concrete algorithm that can enable us to compute Bessel functions with purely imaginary order 
as they have had no applications in areas of natural sciences till date. Only a hand full of pure mathematicians made 
attempts to investigate such functions. 

The author of one of the most detailed treatment on Bessel functions pj thought that these functions were of no 
interest though they were noticed by Lommel ||j who defined the function J v +ifi(%) in the form 



J ^ {x) = r(, + z M + i/2)r(i/2) \ K »,M + iS v ,M\ (5) 

with K and S being real valued functions. Lommel was motivated by the following differential equation: 

x 2 y" + (2a - 2/3u + l)xy' + [a(a - 2/3v) + f3 2 j 2 x 2(3 ] y = 0. (6) 
If any real j3 ^ then equation (|6|) has a solution 

V = x^~ a [AMyx*) + BJ-^x )] . (7) 
Comparing equation ([7]) with the equation 

x 2 y" + axy' + (b + cx 2f3 )y = 0, (8) 

Lommel found a relationship between the coefficients of equation |0) and those of equation (|9|) and expressed them 
in the form: 



/ 9i/-a = -(o-l)/2, (9) 

pv = ^/[(a-l)/2] 2 -fe, h = sTc- (10) 

It turned out that the real valued coefficients (a, b, c) produced in a general case solution of equation ^ with a 
complex index v. At the end Lommel obtained a very cumbersome solution of equation ([5]) as multiples of a complex 
valued functions. He did not propose a method for the calculation of the real valued functions K and S from ([5]). 
Another mathematician Bocher encountered the modified Bessel function with pure imaginary index while he ponder- 
ered over the solutions of Laplace's equation using the method of variable separation in cylindrical coordinate system 
Q. As opposed to Lommel [J, Bocher proposed real valued solutions with the aid of real valued series without 
studying the convergence of these series. Each of the term of Bocher's series was found to be a ratio of polynomial 
functions in is, with coefficients that increased n the order nl.This possed a major problem for the calculation and 
computation via his representation. 

Finally the McDonald's function Ki T (x) also having pure imaginary index was extensively used in its integral form 
to obtain solutions of the Laplace and wave equations for different boundary value conditions by the Soviet mathe- 
maticians M. Kantorovitch and N. Lebedev 0]. They used the integral representation 



oo 

Ki T (x) = J cxp (— x cosh t) cos (rt)dt, x > 0. (11) 



This representation are convenient only for sufficiently large values of the argument x. Unfortunately, for small values 
of x this integral representation is not quite effective, and this is the case where Bessel functions with pure imaginary 
indices become applicable and useful in charge particle dynamics [2| and strain field investigations in nanostructures 

a. 

The above summary shows that there is no effective algorithm for the computation of Bessel functions with pure 
imaginary indices. In the next sections, we provide this algorithm. 



II. METHOD AND CALCULATIONS 



So, Bessel functions with pure imaginary indices are solutions of the equation 

x 2 y" + xy' + {x 2 + v 2 )y = 0, (12) 

where v is a real number. The pure mathematician George Boole was not interested in Bessel functions but was the 
first person who, over a century and a half ago developed an operational method for solving differential equations. 
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Nevertheless in a small paragraph of his article Q , equation (TT2]) was mentioned and the substitution that facilitated 
the solving of the problem of calculating Bcsscl functions with purely imaginary order was used. So expression @ 
suggests the form of the solutions of equation (|12|) : 

y = A (x) cos (y In x) + B (x) sin [y In x) , (13) 
where A (x) and B (x) are series with the coefficients as shown: 



A (x) = a + a\X + a 2 x 2 + • • • + a n x n + ■ ■ ■, 

B (x) = b + b lX + b 2 x 2 + ■■■ + b n x n + ■■■. U4j 

Substituting equation (fT3|l into equation (fT2|) and equating the coefficients before the linearly independent functions 
cos (vhix) and sin (yhix) to zeros, firstly we obtain the relation 

a x = &! = (), (15) 

and secondly Vn, > 2 the recurrent relationships 



na n - 2 - 2^6„_ 2 

( 2 i a 2 \ ' ( 16 ) 

n (n^ + 4^ z ) 

-2va n - 2 + nb n -2 



n (n 2 + Av 2 ) 



(17) 



It is clearly seen from equations (fl7o|) and (jTTJ) that odd term coefficients take on zero values while those with even 
terms can be obtained via a recurrent relations using the values ao, 6o. The solution of the Bessel equatio n |T2|) in the 
form (|13p was firstly introduced by Boole in Q. However, he did not analyze the convergence of the series (fT4|) . Boole's 
recurrent relationship was completely forgotten although they can be easily simplified and used to effectively calculate 
the functions under examination. Since only coefficients with even indices occur, we can re-define the functions A (x) 
and B (x) in the form: 



^h]>> 2 „2 2 ™(!) =E^«(f) ■ ( 18 ) 

B0r) = £&2„2 2 "(f) =£**•(!) • (19) 



n=0 n=0 

Now for n > 1 the recurrent relationship for the coefficients Am , B2„ can be simplified: 



A 2n = ' 2 " - , (20) 

vA 2 n~2+nB 2 n-2 , n1 v 
i>2« ■ (21 J 



III. CONSTRUCTION OF COMPUTABLE SOLUTIONS AND PROOF OF CONVERGENCE 



In order to investigate and numerically compute the strain field distribution in nanostructurcs, or to numerically 
implement these Bcsscl functions in other areas of physics, we need to be able to demonstrate that the solutions to 
this differential equation exist and above of all converges in our domain of definition and interest. We now show that 
the series (fTOjl . ([T7]) converges and does so absolutely for any complex valued argument x and order v. As a series 
will converge absolutely if and only if the absolute value of the n th term (which we shall refer to as the majorant) 
converges. Let the majorant M 2n be defined such that 



M 2n = \A 2n \ + \B 2n \ 

It follows from (HJ) , $TZj) that 



(22) 
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M 2n < 



n 2 + v 2 n (n 2 + v 2 ) 



M 2r , 



< 



M 2 „_ 2 . (23) 



The last inequality provides us the possibility of studying the order of which the majorant M 2n decays (or grows). It 
can be easily shown that 



M 2n < C- 



(n!) 2 



(24) 



where C is a positive constant dependent on the zero-th terms of the sequence. Using the d'Alembert's test for 
convergence the proof of absolute convergence for all values of x and v, is completed. 

Now it is possible to form, without lost of generality, two linearly independent solutions of (fl"2"|) . Let us choose two 
pairs of the values for Aq,Bq that generates the sequences A 2n ,B 2n and its corresponding functions (fl8|) and (|19|) . 
For brevity, these two pairs provide two solutions (|12p that can be easily computed numerically for a few number of 
iterations. In [2], it was shown that the solution spawned from the pair 



(A), Bo) = (0,1), 
represented as Sf„(x) and the ones that are spawned from the pair 



(25) 



(A ,B ) = (1,0), 

was also represented as Cf v (x). Then an analytical expression of these functions can be easily obtained: 



(26) 



Sf„(a;) 



1 - 



1 



1 + v 2 V 2 



+ 



sin (v\tlx) 



1 + v 2 V2 



cos(z/lnx), (27) 



Cf„(a:) 



1 + v 1 V 2 



cos (yhix) 



For x — > we have two equivalent relationships 



1 + v 2 V 2 



sin In a:). (28) 



Sf„(x) ~ sin (yh\x) , 
Cf „ (x) ~ cos (y In x) 



(29) 
(30) 



Now the general solution of the differential equation ([12]) with real coefficients may be explicitly written out via the 
real valued functions 



y = dSf v (x) + c 2 Ct v {x) , (31) 

where c\ and c 2 are any real constants. Provided that equation (|12p has a complex valued solution in the form of 
the Bessel function (|4]). Therefore the latter must be the linear combination of the functions Sf„(x) and Ci v (x). It is 
easy to show that 



Jiv{x) = Cf„(x) + iSfy(x) . 



(32) 
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In other words, the functions Sf„(x) and Cf„(x) are the real and imaginary parts of Bcssel function ([4]). It follows 
from the following equivalent relationships: 



Jiu{x) 



T{iv + 1) V2 



1 - 



1 /x\ 2 



iv + 1 V 2 



~ . v [cos(z^lnx) + i sin(i/ In a;)] ~ [cos(flnx) + isin(i/lnx)] . 
The first equivalent relationship is valid as long as x — > and the second^ one is valid so long as v — > 0. 



As 



V' + ^' + (-x 2 + i/ 2 )y = 0. 



(33) 



cos (i/ mix) = cosh cos { v lux) — sinh (~2~) sm lux) 



(34) 



and as 



sin (y In ix) = cosh ( — J sin (y In x) + sinh ^ 



— i cos (i/lnx) . 



(35) 



C{x) = A{ix) = Y d C 2 J~) 



D(x) = B(ix) = £z) 2 „(f 



n=0 

it is possible to construct two linearly independent solutions for this case in the form 

yi = C (x) cos (y In x) + D (x) sin (z/ In x) , 
j/2 = D (x) cos (i/ In x) — C (x) sin {y In x) , 

where as earlier mentioned, for n > 1 the recurrent equations for C 2ri ,T) 2ri takes the expected form 

nC-in-i - vD 2n -2 



C 2n — 

D 2n = 



n (n 2 + is 2 ) ' 
vC 2n -2 + nD 2n -2 
n (n 2 + ^ 2 ) 



(36) 
(37) 



(38) 
(39) 



(40) 
(41) 



As the series C(x) and D{x) have the same majorant (|24|) , these series will thus converge absolutely for any x and v. 
One real valued solution of (|12|) spawns two real valued solutions for (|33|) . Thus by cho osing only one pair of values 
Co, Do it is possible to obtain two real valued functions being solutions of equation (|33|) . Let 



(Co, A)) = (0,1) 



(42) 



Substituting (f36|) and (f37f into |38|) and ([39]) it is possible to obtain two functions which was denoted in [2| as Sdj, (x) 
and Cd„ (x): 



Sd y (x) 



— (-)' 
l + v 2 \2J 



sin (i/lnx) 



/i\ 2 
^ V2, 



cos(j/lnx), (43) 
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Cd, (a:) 



1 fx\ 2 



1 + I/ 2 V 2 



cos (z^lnx) 



1 + v l V 2 



As x — > wc also have two equivalent relationships 



It is obvious that 



where x and v are any real numbers. 



Cd y (x) ~ cos (v In x) , 
Sdjy (x) ~ sin (i; In a;) . 



Ji V (ix) = Cd v (x) + iSdv (x) 



sin (y In 2) . (44) 



(45) 
(46) 



(47) 



IV. WRONSKIANS OF FUNCTIONS 

It can be easily shown that the Wronskian W of the two linearly independent solutions of the equation 

y" + -y' + f(x)y = 

x 



(48) 



has the form W = a/x, where the constant a depends on the pair of concrete solutions. From the equivalent 
relationships (23]), O and (gSJ), |@BJ) it follows that 



Cf„ (x) [Sf„ (x)]' - Sf„ (x) [Cf„ (x)]' = -, 

x 



Cd„ (x) [Sd„ (x)]' - Sd„ (x) [Cd„ (x)]' 



It is clearly seen that the Wronskians (|49|) and (|50|) vanishes for v = 0. As f — > we have 



(49) 
(50) 



and 



Cf v (x) -> J (x) , 
Sf„ (x) 



Cd^ (x) -> Z (a:) , 
Sd„ (x) -► . 



(51) 
(52) 



(53) 
(54) 



V. ACCURACY OF PROPOSED COMPUTATIONAL ALGORITHM 

To obtain a measure of the computational error for the proposed procedure, we require a detailed study of the majorant 
(remainder) M<x n . Let 



M 2n = ——?m 2 r. 



(niy 

Now it is possible to obtain an estimate of remainder m 2n for n > 2 



(55) 



m 2n 

















n J 





(56) 



The Lagrange's formula for remainder of a Taylor's series gives for n > 2 



m2n <x _ v _l HGH- 1 ) , 
6n 3 



™2n-2 

x 3|i/| + (M-2)(l + M 



1 - 



kl-3' 



where < 8 < 1. The last inequality may be improved by removing from the right hand side the quantity 9: 



TOOn i> F 

— < 1 1 

TO2n-2 n 2 n 3 



where 



^IH-l^-M) |i/|<2, 
F = <J |H |H - 1| [3\u\ + (\u\ - 2)(1 + M/2)2 3 -M] /6, 2< |i/|<3. 

||z/| - 1| {v 2 /2 + Z\v\ -2)/6, > 3. 



As for v ^ 



then it is possible to obtain 



rn-2-n v 2 F v 2 F 

— — < 1 1 < H 1 , 

rri2n-2 Ti a n z n a 



v 2 F v 2 F 
In (m 2n ) - In (m 2 „_ 2 ) < ln(l + + -^) < + -s • 

n n n n° 



Using the well known sums 



00 -i 2 

E J 2=T- 1 - ' 6449 ' 

n— 2 

oo - 

£ -j = C(3) - 1 = 0.2021 



n=2 



and also the values 



m = 1 , 

i+M 



TO 2 

it is possible finally to obtain the inequality 

In m2n < 0.6449^ 2 + 0.2021F . 

TO2n-2 

From the last inequality it is easy to obtain a final one simplified in the form 



m 2n < ] + \ V \ exp(0.6449w 2 + 0.2021F) = m(v) . 

1 + IS Z 



Now we have 
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M 2n <m(i/)— 2 
(n!) 



(63) 



and it is possible to conclude that the sum of the remaining terms for the series of functions A (x) , B (x) , C (x) , D (x) 
cannot exceed the value 



Since the remainder depends on the value of v, let us consider for example the case \v\ < 2. Then 



(64) 



£at < m 



x 2 °° 1 



2 00 1 2 

4)jEifT- (65) 



-N 



(n!)"V2 



The sum at the right hand side of equation (|65|) is simply the remainder of the series for the modified Bessel function 
Jo (x). This remainder may be evaluated further as shows below: 



1 fx\ 2n 1 /x\ 2N 



£ (»» a (2) 



/xy jv ^ (-W!) (x\ 2n 
V2J ^ [(n + iV)!] 2 

1 /x\ 2N 



(my 



(my 



1 



(l + JV)" 



[(l + iV)(2 + A0--.(n + A0]' 



(!) 



2// 



< 



< 



1 /X\ 2N 



Now for any real x and M < 2 we have 



1 /x\ 2 



1 • 2 V2 



n! (n + 1) 



I ( X -\ 

+ 1)! 12/ 



1 /x\ 2jV 2 
(TV!) 2 v2/ 



1 z^^- 1 



Ji (a?) • (66) 



I 1 Iwl 1 /x\ 2iv +l 

< 1 ' cxp (0.6449i^ 2 + 0.2021F) ? ( - ) /j (x) 

~ 1 + z/ 2 v ; (iV!) 2 V2/ 



(67) 



Elimination of |z/| in the last inequality gives an estimation for any x and \v\ < 2: 



^^(2) ^ 

For instance for x < 2 inequality (|68|) reduces to 



24 

ejv < o- (68) 

" (AH) 2 

The last inequality shows that for the functions A (x) , B (x) , C (x) , D (x) and for x < 2 and \i/\ < 2 a computational 
error less then e% < 1.5 x 10~ 16 occurs after the computation of only eight terms of the corresponding series. 
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VI. CONCLUSION 



A new algorithm for computing the real valued solutions of some special functions namely the Gamma function with 
pure imaginary argument, the Bcsscl (and modified Bcssel) function of purely imaginary orders was systematically 
declared. We have used these developed Bessel functions for the investigation and calculation of strain fields in 
semiconductor structures @, [H| . In the above applied case the values of the arguments of the Bessel functions lie 
in the vicinity of unity. The algorithm for their computations is shown to be rapidly converging and gives a very 
favorable accuracy for as few as eight iterations. 

The algorithm described may give an alternate approach to solving the ground state of the screened coulombs 
potential problem. It is also promising in solving the inverse scattering problem which is mostly approached via the 
inverse Green's function. Further studies of the asymptotic behavior of these functions are still warranted as they 
play a great role in a wide class of boundary value problems and integral and differential problems of wave scattering. 
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